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We investigate the effect of noise on Random Boolean Networks. Noise is implemented as a prob- 
ability p that a node does not obey its deterministic update rule. We define two order parameters, 
the long-time average of the Hamming distance between a network with and without noise, and the 
average frozenness, which is a measure of the extent to which a node prefers one of the two Boolean 
states. We evaluate both order parameters as function of the noise strength, finding a smooth tran- 
sition from deterministic (p = 0) to fully stochastic (p = 1/2) dynamics for networks with K < 2, 
and a first order transition at p = for K > 2. Most of the results obtained by computer simulation 
are also derived analytically. The average Hamming distance can be evaluated using the annealed 
approximation. In order to obtain the distribution of frozenness as function of the noise strength, 
more sophisticated self-consistent calculations had to be performed. This distribution is a collection 
of delta peaks for K = 1 , and it has a fractal sructure for K > 1 , approaching a continuous distribution 
in the limit K > 1 . 

PACS numbers: 89.75.Da,05.65.+b,91.30.Dk,91.30.Px 
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I. INTRODUCTION 

Random Boolean Networks (RBNs) 0] have been 
used as a simple model for a variety of dynamical sys- 
tems consisting of interacting units, such as neural net- 
works 01, social networks 0] and, more prominently, gene 
regulatory networks 00]. RBNs are composed of Boolean 
nodes that are coupled to each other. In the case of gene 
regulatory networks, the Boolean state is a step-function 
approach to the expression level of a particular gene. De- 
spite this loss of detail, the most important features of gene 
regulatory processes are still captured in many cases, since 
they should not depend on biochemical details, but on the 
desired sequence of events in the cell 0]. 

So far, the dynamics of RBNs have mostly been studied 
using deterministic update rules. The dynamics of such 
models are non-ergodic, with periodic attractor trajecto- 
ries in state space. Once the system has reached an attrac- 
tor, it remains there. Another important property of RBNs 
is a phase transition, which occurs when the number K of 
inputs per node is changed. For unbiased networks, the 
dynamics exhibit a frozen phase at K = 1 , where local per- 
turbations die out quickly and most attractors are fixed 
points, and a "chaotic" phase at K > 2, where perturba- 
tions increase exponentially fast and attractors have very 
long periods. At the boundary K = 2 between those two 
phases are the so-called "critical" networks, where pertur- 
bations increase algebraically with time. Originally, it was 
suggested by Kauffmann 0] that such critical networks are 
best suited to model real systems, which are supposedly 
poised "at the edge of chaos". In the meantime, there is 
agreement that RBNs of all three types have only limited 
validity when applied to real systems. 

Real networks usually have some level of stochastic be- 
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haviour, and for this reason several authors have investi- 
gated RBNs under the influence of stochasticity. For in- 
stance, in 0] the nodes of RBNs were updated in a com- 
pletely random order. This update method preserves the 
non-ergodicity of the system, and it is still possible to 
identify distinct attractors. Attractors are in this case de- 
fined as sets of states all of which are visited for a non- 
vanishing proportion of time during the same trajectory. 
The stochastic update sequence vastly reduces the num- 
ber of attractors of critical RBNs, which becomes a power 
law as function of the network size. Similar results are 
obtained when the update sequence deviates only slightly 
from a synchronous update |cS|. Such a power law was for 
a long time falsely believed to occur in deterministic RBNs 

Sim El. 

Instead of introducing stochasticity into the update 
times, other authors introduce it into the update functions. 
In (H|], probabilistic Boolean functions are used, where a 
set of several Boolean functions is assigned to each node, 
and at each time step one of these is chosen randomly with 
a given probability. According to fl^l . this model is more 
realistic than models with a purely deterministic update 
scheme. 

However, the most important way of introducing noise 
into a RBN is in form of a "temperature", leading often 
to ergodic behavior. The effect of thermal noise on Ising 
spins on a network was studied in Il3l Hill , where a "fer- 
romagnetic" transition from the ordered to the disordered 
phases was observed at a critical noise strength value. In 
the language of gene regulatory networks, a temperature 
manifests itself as fluctuations in the protein concentra- 
tions, so that a gene may not always be turned on or off, 
given the same expression state of the other genes fl5ll . 
This effect can be included into models by allowing a de- 
viation from the deterministic update rule with a certain 
probability. In jloTl , for instance, a subset of nodes were 
perturbed in this way (this corresponds to turning on the 
temperature for a short time interval), and the response 
of the dynamics to this perturbation was evaluated, giv- 
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ing information about the basin structure of the system. 
Miranda et al. fl7ll studied the effect of a permanently act- 
ing temperature by introducing a fixed probability p that 
the state of a node becomes the opposite of what it should 
be according to the deterministic update rule. They evalu- 
ated the average crossing time between trajectories in state 
space which started from different initial states as function 
of noise strength il^,ITR[l9ll . By sampling the entire state 
space of small networks (N ^ 20), it was found that the 
"barriers", which correspond to the attractor basin bound- 
aries, can be crossed with non-vanishing probability when 
p > 0, although the characteristic times may be large. This 
means that the system is always ergodic. This type of noise 
has also been studied for Boolean networks with threshold 
functions, corresponding to a majority update rule l2(3l . 
This system undergoes a second order phase transition at 
a critical noise strength from an ordered dynamical phase, 
where all nodes assume the same value for the majority 
of time, to a disordered phase where nodes assume both 
states equally often. 



In this work, we investigate the effect of ongoing 
stochastic noise on RBNs. Following JUL M, B, noise 
strength is tuned via a probability p that a node does not 
obey its deterministic update rule. We monitor the transi- 
tion from fully deterministic dynamics (p = 0) to purely 
stochastic dynamics (p = 1 /2) as the noise strength is var- 
ied. Differently from fT^,ll1lll9tl , we are interested in the 
behaviour of the networks in the limit of large system size, 
where it is impossible to explore large parts of the state 
space. In order to characterize the transition from zero 
to infinite temperature, we define two order parameters: 
the long-time average of the Hamming distance between a 
network with and without noise, and the average frozen- 
ness, which is a measure of the extent to which a node 
prefers one of the two Boolean states. We find, both ana- 
lytically and numerically, that this transition for the Ham- 
ming distance is continuous for K ^ 2, and discontinuous 
at p = for K > 2, when the Hamming distance is consid- 
ered. This distinction is a direct consequence of the phase 
transition from frozen to chaotic dynamics in the deter- 
ministic model. The frozenness shows a smooth transition 
for all values of K. The distribution of frozenness shows a 
surprising richness in structure, as revealed by computer 
simulations. For K ^ 2 and for K » 1, we succeeded in 
reproducing this structure as function of p by analytical 
considerations. 



The remainder of this paper is divided into the follow- 
ing parts: In Sec.|n]we define the RBN model and the type 
of noise used for our study In Sec. [TTTJ. we define the first 
order parameter, the Hamming distance, and evaluate it 
numerically and analytically. In Sec. lIVl we define the sec- 
ond order parameter, the frozenness, and evaluate it using 
computer simulations and analytical considerations. Fi- 
nally, we summarize and discuss our findings in Sec.lVl 



II. MODEL 

A Boolean network is defined as a directed network of 
N nodes representing Boolean variables er e { 1 , 0} N , which 
are subject to a dynamical update rule, 

<r(t + 1) = /(«r(t)), (1) 

where f i is a function assigned to node i that depends ex- 
clusively on the states of its inputs. 

We introduce noise into the system through a probabil- 
ity p that a node does not obey its deterministic update 
rule, 

<x(t + 1) =/(<x(t))y n , (2) 

where n is a random vector, with elements being 1 with 
probability p and otherwise. The symbol Y. represents 
the "exclusive or" Boolean operation. Hence, for p = the 
deterministic behaviour is recovered, and for p = 1 /2 the 
dynamics is completely stochastic. 

RBNs are a special case of Boolean networks, where all 
possible Boolean functions are assigned randomly to each 
node with the same probability, and where the nodes are 
randomly connected. The number of inputs of each node 
is fixed at a value K. The random wiring leads to a Pois- 
son distribution with mean K for the number of outputs. 
When updated deterministically RBNs are in the frozen 
phase for K = 1 . After a transient time, they reach an at- 
tractor where all nodes (or all nodes apart from a small 
number) are permanently frozen in one of the two Boolean 
states. Networks with larger K have also a frozen core of 
nodes for p = 0, and the nodes belonging to it become 
frozen after a transient time. For K = 2, all but of the order 
of N 2 / 3 nodes belong to the frozen core. With increasing 
K, the frozen core contains an ever smaller proportion of 
nodes. For K = 2, the nonfrozen part of the network con- 
sists of several independent components. Each of these 
components contains a set of relevant nodes, which are con- 
nected such that there is at least one feedback loop among 
them, and "trees" of nonfrozen nodes which are rooted in 
the relevant nodes and which are slaved to the dynamics 
of the relevant nodes. 



III. AVERAGE HAMMING DISTANCE 
A. Definition 

We use the average in time of the Hamming distance 
between the states of two copies of a network in order to 
quantify the effect of noise on the dynamics. Consider a 
given network in the initial state <r(t = 0), and an exact 
replica, which is initially in the same state, <r'[t = 0) = 
er(t = 0). The dynamics of both networks are evolved in 
parallel, but noise is applied only to er'(t), as in Eq. ((2). 
The mean Hamming distance h(t) between the two net- 
works is defined as 

h(t) = l(|ff l (t)-a((t)|> / (3) 
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where (...) denotes the average over the noise. 

The long-time average h of the Hamming distance is de- 
fined as 

h= lim 11 V | ai (t)-a((t)|. (4) 

T— >oo I IN * 

t/l 

If the trajectories become completely uncorrelated after 
some time, we have h. = 1/2. If the trajectories remain 
closer in state space, we have h < 1/2. The case h > 1/2 
does not occur in our model and is therefore not consid- 
ered in this paper. 



several values of K. The solid lines are the fixed point so- 
lutions of Eq. ©, the symbols represent the result of com- 
puter simulations of quenched RBNs. The agreement be- 
tween the annealed approximation and the real networks 
is very good. 

The most striking feature of Fig. [TJ is the existence of 
a first-order transition at p = for K > 2. This is due 
to the phase transition to "chaotic" behaviour for K > 2. 
In chaotic networks, even the smallest local perturbations 
have a global effect. 



C. The Hamming distance on subsets of nodes 



B. Annealed approximation 

We will first evaluate analytically the Hamming dis- 
tance by using the so-called annealed approximation |2lll . 
This is a mean field theory, which neglects correlations 
between nodes and the finite size of the network. The 
annealed approximation corresponds to the behaviour of 
a (infinitely large) network where all the edges are ran- 
domly rewired at each time step. Within the annealed ap- 
proximation, the dynamics of a RBN without noise (i.e. for 
p = 0) is fully specified by the parameter A, which is K 
times the probability that a node changes its state when 
one (or more) of its inputs is flipped. For RBNs, we have 
A = K/2, since for any input combination, there is an equal 
probability that the output of a function will be either or 
1 . When considering two replicas of a network, A is iden- 
tical to the mean number of nodes that assume a different 
state in the two networks at time t = 1 when at time t = 
the state of only one node was different. 

At any time, the Hamming distance between a network 
with noise and its twin noiseless counterpart, as described 
by Eq. (0, is simply the fraction of nodes which were 
changed by noise or by the effect of previously changed 
nodes. The time evolution of h(t) can then be described as 
the evolution of the population of flipped nodes. (A node 
in the replica with noise is called "flipped" it its state de- 
viates from the state it has in the replica without noise.) 
Let q(H(t)] = 1 - (1 - h(t)) K denote the probability that a 
node has at least one flipped input. Then the probability 
H(t + 1 ) that a node is flipped at time t + 1 can be written 
as 



We next evaluate separately the Hamming distance for 
the frozen core and for the nonfrozen part of the network. 
Fig. [2] shows the long-time Hamming distance, evaluated 
only for the nodes that belong to the frozen core. These 
curves can be fitted using the annealed approximation 
Eq. 10 under the condition that the factor A/K on the right- 
hand side of Eq. I0 (representing the probability that a 
node is flipped when at least one input is flipped) is re- 
placed with A e ff/K, with A e ff being used as a fit parameter. 
For K = 1 and 2, the frozen core is virtually indistinguish- 
able from the rest of the network, and A e ff = A, but for 
K > 2, A e ff decreases with increasing K. The reason is that 
the frozen core becomes composed mainly of nodes with 
constant functions which have A e ff = 0. 




h(t + r 



A 
K 
AO 



q(1-p] 
2p) 



pq + (1 -q)p 



(5) 



K 



l-O-Ti(t)) 1 



+ V, 



where the first term in the first line corresponds to the pro- 
portion of nodes that are flipped by previously flipped 
nodes (and are not flipped back by noise), and the sec- 
ond and third term are the proportion of nodes that are 
flipped by noise (with or without inputs being flipped). 
The fixed point of Eq. 10 determines the order parameter 
h, for given K and p. We evaluated this fixed point numer- 
ically. Fig.[T]shows h as function of the noise strength p for 



FIG. 1: (Color online) Average Hamming distance as function 
of noise strength for RBNs of size N = 1 4 for different values 
of K. Each point was obtained by averaging the results over 3 
different network realizations. The solid lines are the respective 
steady-state solutions of Eq. $5$. 

Before evaluating h for the nonfrozen nodes, let us con- 
sider the simplest possible connected set of nonfrozen 
nodes, which is a simple loop. For K = 1 and K = 2, 
such simple loops of nonfrozen nodes play an important 
role at determining the attractors with periods larger than 
1 l22ll , however, a considerable fraction of K = 2 networks 
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also have more complex relevant components. The effect 
of noise on such loops is very different from its effect on 
the frozen core, since if one of its nodes is flipped, this flip 
propagates indefinitely around the loop. One can calculate 
the accumulation of flips on such loops by considering the 
average Hamming distance at a given time between a loop 
without and with noise, 



h(t)= y_ 

i=0 



i=0 



21+ 1 

(2i+1)! 



P 2i+I d-P 



it-(2i+1) 



(6) 



1 



1 - e- 2t ? 



The first equation evaluates the probability that a node has 
been flipped an odd number of times, and the subsequent 
transformations are valid for t > 1. The Hamming dis- 
tance approaches the value 1 /2 with an exponential decay, 
and with an characteristic time T = 1 /(2p). 




K = 1 , A e ff = 0.50 
K = 2, A eff = 1 .00 
K = 3, A ef f = 0.33 
K = 4, A ef f = 0.00 



K 2, the data points are considerably above this expo- 
nential curve because a node can become flipped via many 
different paths. The data are better fitted using Eq. (O, in 
particular for long periods (i.e. large times). Just as for the 
case of the frozen core, A e ff was used as a fit parameter. 

For smaller values of the attractor period, the data are 
considerably below the fitted line. The reason is that these 
attractor periods are much smaller than typical attractor 
periods, and networks with such short attractors are not 
characteristic of the ensemble, but have a state-space struc- 
ture with a smaller set of recurrent states. Consequently 
trajectories diverge less fast than in typical networks. 




FIG. 3: (Color online) Average Hamming distance h(l) of the rel- 
evant components of RBNs with different values of K, after a full 
period I, with p = 0.01 . The curves were obtained by sampling at 
least 2 x 1 4 attractors of several distinct RBNs of sizes N = 1 2 , 
50 and 25 for K < 2, 3 and 4, respectively. The solid lines are 
given by Eq.[6]for K > 1 and Eq.|5]for K > 2. 
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FIG. 2: (Color online) Average Hamming distance as function 
of noise for RBN of size N = 1 4 for different values of K. Only 
nodes belonging to the frozen core of the network (without noise) 
were considered. Each point on each curve was obtained by av- 
eraging the results for 3 different network realizations. The solid 
lines are the respective solutions of Eq. {5), with A ef f being used 
as a fit parameter. 

We evaluated how fast a trajectory leaves an attractor in 
the presence of noise by first letting the system approach 
an attractor and by then turning on the noise and measur- 
ing the Hamming distance h(l) to the initial state after one 
attractor period I. This is identical to the distance from the 
state of the noiseless replica, which returns to the initial 
state at time I. Fig. [3] shows the values of h(l) for RBNs 
with different values of K. For K = 1, the data match 
Eq. ((6) very well, since the nonfrozen part of the network 
in this case can only be composed of simple loops. For 



IV. FROZENNESS 
A. Definition 

The "frozenness" of a network measures the extent to 
which the nodes spend more time in one of the two 
Boolean states. It is zero, when the nodes spend the same 
time in both states, and it is 1 when the network is frozen. 
The frozenness of node i is defined by the expression 

a i = (qW- q W) 2 , (7) 

where q& ^ is the proportion of time node i is in state a, 

1 T 

q&° = Km =Y- b °dt).<r ( g ) 
1 t=o 
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By eliminating one of the two probabilities from Eq. (0, 
we obtain 

Qi = (iq^-tf, (9) 

where cr is either or 1 . 

The frozenness (CI) of the network is obtained by aver- 
aging over the nodes. Figure H] shows the frozenness (CI) 
as function of the noise strength p obtained by computer 
simulations of networks of size 10 4 , for different K. 
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FIG. 4: (Color online) Frozenness as function of noise strength 
for RBNs of size N = 10 4 for different values of K. Each 
curve was obtained by averaging the results for 100 different net- 
work realizations. The solid lines correspond to the averages of 
p n (0|K,p), obtained in Sees. HVBl HVCl and HVDl for K = 1, 2 
and > 3, respectively. 

At p =0, the frozenness corresponds obviously to the 
size of the frozen core, and it decreases towards as the 
noise strength approaches the value 0.5. 

In order to derive these curves analytically, the annealed 
approximation is of no use, since a network that is rewired 
during the course of time has very small frozenness, which 
is due uniquely to the constant functions. Therefore, a sim- 
ple analytical calculation, which does not require the con- 
sideration of correlations between nodes, can only be per- 
formed for nodes with constant functions: for such nodes 
the frozenness is given by 

Qi-(l-2p) 2 . (10) 

In the following, we will present more advanced ana- 
lytical evaluations and further computer simulations for 
RBNs with different values of K. 



B. K = 1 

For K = 1 there are only 2 2 ' =4 possible Boolean func- 
tions, two of which are constant (1 or 0), and the remain- 
ing ones are the copy (f(o~) = a) and invert (f(cr) = ^o~) 



functions. As far as the analysis of frozenness is con- 
cerned, there are only two distinct functions, constant and 
non-constant, since the output value is not relevant, but 
only how often it changes. Since each of the two types of 
functions occurs equally often in a K = 1 RBN, we have 
A = 1 /2, but K = 1 networks with other values of A can 
also be constructed. 

In a network with K = 1, each node has only one in- 
put, and this input node also has one input, etc. In order 
to evaluate the probability that a node is flipped, one only 
needs to consider the chain of those nodes that can have 
an influence on the considered node. Nodes with constant 
functions present a barrier to the propagation of a pertur- 
bation, since they do not respond to a change in their in- 
puts, and therefore the chain ends (or, more precisely, be- 
gins) at a node with a constant function. 

Without loss of generality, we define q' 1 ' as being the 
proportion of time node i assumes its most frequent value, 

Emaxfq^l-q? 1 ) e [1/2,1]. (11) 

Since the value of q' 1 ' is fully determined by the distance 
of node i to a node with a constant function, we choose 
the label i in the remainder of this subsection to signify 
this distance. If the node itself has a constant function, we 
have i = 0, if the node has a non-constant function, but its 
input has a constant function, we have i = 1, etc. 

The value of q ' ' , i.e. for nodes with constant functions, 
is simply 

q<°>(p) = 1-p. (12) 
For larger values of i, we have the recursion relation 
qW = (1-p) q H) + (i_ q (i-l)) p 



where the solution of the recursion relation was obtained 
using Eq. Q2). The probability of finding a given q' 1 ' in 
the network is 

Pq(q (i) |p) = (l-A)A i = Qy + . (14) 

The frozenness of the network is thus given by 

(a)[ V )=Y_ Vci [q^\ V )^-2q^) 2 / (15) 

i 

which is plotted in Figure |4] and fits the curve for K = 1 
well. 

Although networks with K = 1 have only a discrete 
set of possible q-values, the distribution of q values ap- 
pears as a continuum when determined by computer sim- 
ulations. There are two reasons for this. First, the data 
points for q close to 1/2 (i.e. for CI close to 0) are so 
close to each other that they cannot be resolved, since a 
computer simulation uses a non-vanishing bin size. Solv- 
ing Eq. (TT3T) for i, inserting the result in {141 , and using 
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the relation p q (q|K = 1,p)dq = p q (q|K = 1,p) with 



dq 



, we obtain 



p q (q|K = 1,p) oc (2q-1)tad-2p) 



(16) 



with A = 1/2 for RBNs. Thus, in the limit q -> 1 /2 the dis- 
tribution of 2q — 1 follows as a power-law with an expo- 
nent given by the above expression. The probability den- 
sity of D. also decays as a power law, in the limit Q — > 0, 
but with a different exponent, since 



Pa(H|K = l,p) = p q 



K = 1,p ^ 

W da 



(17) 



QIn(1-2p) 



-3/2 



Second, a computer simulation averages only over a fi- 
nite amount of time, T, and therefore the measured values 
q ' are Gaussian distributed around the exact value q, 



p q ,(q'|K = l,T,q) 



1 



(18) 



This dependence on T can be included in Eq. (l4|l to ob- 
tain the probability density function for q 



p q (q|K = 1,p,T) =(1-A)^^Pq'(llT,q (i) 



i=0 



For the distribution of £1 values, we obtain 

'Vn + i 



Pa(0|K,p,T) = p p 



K,p,T 



2/Q' 



(19) 



(20) 



Fig. [5] shows the distribution of the frozenness for a 
quenched network with K = 1, for several values of p. 
It can be seen that there is very good agreement with 
Eq. | |20l l. The presence of fluctuations significantly devi- 
ates some of the distributions from the expected power- 
law decay. For p = 0.01, the small values of frozenness, 
which are not in agreement with the theoretical result, 
are due to the existence of loops, which are omitted in 
the analysis above. The probability that a node is part of 
a nonfrozen loop tends to zero as the network becomes 
larger, and therefore these points vanish in the limit of in- 
finite system size. 



C. K = 2 

As K becomes larger, the number of possible functions 

grows very fast as 2 2K , and a detailed analysis of the 
frozenness, as was done for K = 1, becomes more com- 
plicated. The values of q are still discontinuously dis- 
tributed, but their number increases fast with K, due to 
the numerous combinations of Boolean functions that can 
determine the q -values of the K inputs of a node and thus, 




FIG. 5: (Color online) Distribution of the frozenness for differ- 
ent values of the noise strength for RBNs of size N = 10 5 for 
K = 1 and T = 10 4 . Each curve was obtained by averaging over 
100 different network realizations. The solid lines are given by 
Eq.{20}. 



in combination with the node's Boolean function, the q- 
value of its output. Here we will lay out the basic con- 
siderations needed to obtain the distribution of q for all 
K > 1 , and we will obtain by numerical iteration the dis- 
tribution for K = 2. Without loss of generality, we redefine 
q as q = q CT= i , i.e., the fraction of time a given node has 
the value 1 (as opposed to Eq. (jll) . which simplified the 
case K = 1). 

In general, the probability density function p q (q|K = 
l,p) needs to account for all possible recursive combina- 
tions of output functions and their inputs. We can thus 
write the following self-consistent expression, 



p q (q|K,p) 



2> f &(q< f '-q)np q (q (i) i K 'P) d q (i) ' 



1=1 



(21) 



r 



where the sum is taken over all Boolean functions; Pf is and 
the probability of the fth Boolean function (pf = 2 _z ), 



.(f) - 



l-2p)q( f }({q (i) })+P, 



(22) 
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TABLE I: Expressions ofq' f '(qi,q2) for all Boolean functions for 
K = 2. The Boolean expressions of each function is also given for 
reference. 

where q ' f ' ({q ' x '}) is the value of q for a specific function 
f , given the values {q } of its inputs, for i = 1 , . . . , K. 

Since Eq. d211 involves an expression q' f ' ({q' 1 '}) for all 
Boolean functions, a general closed solution becomes un- 
feasible. However, for K = 2 Eq. lT2T) l can at least be 
solved numerically, since there are only 1 6 possible func- 
tions, given in Table |U Eq. iTZTl l is then solved by iteration, 
until convergence to a self-consistent q distribution is ob- 
tained. We started with the initial distribution 

Pq(q|K,p) = ^6(q-p) + l6(q-(1 -p)). (23) 

In the end, we determined the final distribution Pq (Q|K = 
2,p) by using Eq. 420t . 

Fig. [6] shows the distribution of CI for simulated 
quenched RBNs with K = 2 for different values of p, 
compared with the result of the numerical evaluation of 
Eq. jzTlT l as described above. There is a very good agree- 
ment between the two types of results. The peaks corre- 
spond to prominent values of the frozenness. The right- 
most peak is always due to the constant functions, but 
large frozenness values are also obtained for other func- 
tions. For instance, fi assumes the value whenever 
both inputs are different. If both inputs have q = 1/2 
(i.e., Q. = 0), the value of qf = i is (1 — 2p)1/4 + p and 
Q = ((1 - 2p)/2 + 2p - I) 2 . This is the second mam 
peak of Pq (Q|2, 0.4, 1 4 ) (counted from the right end). For 
smaller values of p, the peaks are not discernible, and a 
broad continuum appears, with a distribution that follows 
a power-law with an exponent ~ 0.6 as CI — > 0. As T 
becomes larger, it is expected that the continuous regions 
become more and more discontinuous, as can be seen in 
Fig. which shows the theoretical prediction for larger 



times T. Moreover, when the resolution is increased (see 
inset), it can be seen that peak-like regions which appear 
like fluctuations around a single value of CI, are in fact 
composed of sharper peaks, which themselves are com- 
posed of other peaks, building a fractal structure. 

Another distinguishing feature seen in Fig. [6] is a sharp 
transition at p = 0, where the only two possible values of 
q are and 1 , both of which amount to CI — 1 , leading 
to the variance o"q = 0. For p > this abruptly changes, 
and a wide range of values of q are possible, which discon- 
tinuously leads to o"q > 0. There is no such discontinuous 
transition for other K values, but a continuous one (see fol- 
lowing section), which makes the case K = 2 special. 




Q 

FIG. 6: (Color online) Distribution of the frozenness of nodes for 
different noise strength for RBNs of size N = 1 5 and for K = 2 
and T = 10 4 . Each curve was obtained by averaging over 100 
different network realizations. The solid lines show the values 
of Pq(0|K = 2, p, 10 4 ) according to Eq. l |20l . The line segment 
corresponds to a power law with an exponent 0.6. 

The function (CI) is obtained by performing the integral 
J CI ■ pa [C1\K = 2,p, T)dO. As can be seen in FigurelH our 
calculation of this function agrees well with the results of 
computer simulations. 



D. K > 2 

For larger values of K, numerical solutions of Eq. iTZTT l be- 
come progressively more elaborate. We did not pursue the 
task of writing the expressions of 256 function qf ({q' 1 '}) 
for K = 3 or of 65, 536 functions for K — 4. Instead, we per- 
form in the following an approximation that is good for a 
large number of inputs per node. 

When K is large, the vast majority of Boolean functions 
have the output 1 for approximately half the input com- 
binations. This means that almost all nodes have at their 
inputs q values close to 1/2. We therefore make the as- 
sumption that the input values to each node are 1 and 
with probability 1 /2, independently from each other. This 
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FIG. 7: (Color online) Expected distribution of frozenness of each 
node as function of noise for RBNs with K = 2 and different 
values of T. The curves represent Eq. {20}. The inset shows a 
zoom into pn (0|K = 2, 0.4, T) for different values of T. 



means that for any given function all input combinations 
are equally probable. It then follows immediately that the 
possible q values are identical to the possible fractions of 
output values 1 in the truth table of a Boolean function, 
and that the probability for a given q value is 



p q (q|K,p = 0) = 2 



M 



/ M 
VqM 



(24) 



Here, we have defined M = 2 K , and the possible q values 
are thus multiples of 1 /M. 

In the presence of noise, each output value is inverted 
with probability p, implying that q is changed to q' = 
q(1 — p) +p(1 — q) = (1 — 2p)q +p. We therefore have 



Pq(q|K,p) 



M 



q-p 



K,p = . 



1 -2p KM V 1 -2p 
For the frozenness CI we obtain the distribution 



(25) 



p Q (D|K,p) = 



M 



2(1 -2p)v / n 



Pq 



Vo + 1 -p 



2(1 -2p) 



K,p 



(26) 

Finally, one needs to take into account the effect of fluc- 
tuations, exactly as was done for the previous cases, 



p q (q|p,K,T) = 



p q (q'|p,K)p q (q|T,q')dq', (27) 



where p q (q|T,q') is given by Eq. dl8t . 

Fig. [8] shows the distributions of frozenness for K = 4 
and 5. In contrast to the cases K = 1 and K = 2, the peaks 
are less pronounced, and are hardly visible. For K = 3 (not 
shown) there are some peaks which are still visible, spe- 
cially for high values of p . Therefore, the high-K approx- 
imation Eq $27} is very good already for K = 4. Both dis- 
tributions show the same power-law decay po(fl|K,p) ~~ 



£1 1 / 2 . This is simply due to the fact that for CI — > 
(q — > 1/2) the shape of p q (q|K,p) is essentially flat, and 
thus p n (n|K,p) ~ dq/dH = Q _1 / 2 /2. 




FIG. 8: (Color online) Distribution of the frozenness of nodes as 
function of noise strength for RBN of size N = 1 4 for different 
values of K. Each curve was obtained by averaging the results for 
100 different network realizations. The solid lines correspond to 
Eq.|2Zl 

The quality of our approximation can also be assessed 
by comparing the analytical prediction for (CI) with the 
computer simulations (Fig. [4j. We expect that the approx- 
imation becomes even better for larger K. 



V. CONCLUSION 

We have investigated the effect of thermal noise on 
RBNs by evaluating two order parameters, the long time 
average of the Hamming distance between two networks 
and the average frozenness of the network. While for 
K = 1 and K = 2 the average Hamming distance increases 
continuously from to 1/2 as p increases from to 1/2, it 
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has a jump at p = for K > 2. These findings are well re- 
produced by the annealed approximation, and they are a 
consequence of the transition from a frozen to a "chaotic" 
phase in the deterministic system. In the chaotic phase 
(occuring for K > 2), initially nearby trajectories become 
eventually uncorrelated. The smooth increase of the Ham- 
ming distance towards the value 1/2 is compatible with 
what was found in iT^,IT9tl for small networks. 

The analysis of the average frozenness of the network 
required more sophisticated calculations than the an- 
nealed approximation, and revealed intricate details of 
the network dynamics. For all values of K the probabil- 
ity distribution of frozenness is a set of delta peaks. For 
K = 1, these peaks can be obtained by considering the 
distance of nodes to nodes with constant functions. For 
K > 1, the analysis becomes a lot more elaborate, due to 
the large number of Boolean functions and the resulting 
vast number of possible combinations of frozenness val- 
ues for the inputs of each function. We explained the gen- 
eral method, and performed the actual numerical evalu- 
ation for the case K = 2. The delta peaks show a fractal 
structure, which emerges from the iterated recursion rela- 
tion for the possible frozenness values. The variance of the 
frozenness distribution changes continuously with p for 
all K ^ 2, but for K = 2 it has a jump at p =0, where the 
variance changes discontinuously from to a value larger 
than zero. For larger values of K, the delta peaks are so 
close to each other that the frozenness distribution appears 
continuous, and in this limit we succeeded in performing 
an approximate analytical calculation. 

We do not find a phase transition at finite noise strength, 
in contrast to [20], where Boolean networks with threshold 
functions following a majority rule were used. Such a sys- 
tem undergoes a second order phase transition from an or- 
dered "ferromagnetic" phase, where all nodes assume the 
same value for the majority of time, to a disordered phase, 
where the nodes assume both states equally often. The 
presence of an ordered phase is a direct consequence of 
the majority rule, and this transition is similar to that in a 
network of Ising spins fl3l[T3l . The order parameter in 1 20] 
was defined as the average "alignment" s = |(1 — 2cr)|, 
which is 1 if all nodes are in the same state. The order pa- 



rameter s is only meaningful in systems where the system 
is ordered in the absence of noise, and where the symme- 
try between the states with values and 1 is broken, as 
in ferromagnetic spin systems. Otherwise, (CI) is a bet- 
ter order parameter, because it captures disordered frozen 
phases, such as for K = 1 in RBNs. Of course, a phase tran- 
sition in the value of s is always accompanied by a phase 
transition in the value of (CI) . The opposite is not always 
true. 

It is to be expected that real networks show some kind 
of robustness to noise, since they must be able to carry out 
their function in a noisy environment. As the results of this 
work show, only for RBNs with K = 1 do the order param- 
eters change slowly as noise is switched on. RBNs with 
K > 1 fail to exhibit robustness to noise, which is hardly 
surprising given the random wiring of the system and the 
random choice of functions. It will therefore be interest- 
ing to extend the present study to networks with a more 
restricted set of functions with more biological relevance, 
such as threshold [23] or canalizing functions j23,l25ll . At 
least for some sets of functions, one should expect a phase 
transition at a finite noise strength, similar to the transition 
seen in i20Tl . The survival of the "ordered" phase up to a 
certain noise strength can be viewed as a certain type of 
robustness. 

It remains to be seen how other network topologies l26ll 
and the incorporation of redundancy |2^1 change a net- 
work's response to noise. In j26ll , it was shown that a 
scale-free input distribution changes the average number 
and length of attractors. In |2^1 , redundancy was intro- 
duced as functional duplications of nodes in the network, 
which resulted in greater robustness against random mu- 
tations of the update functions. In both papers, only de- 
terministic dynamics were considered. The effects of these 
(or other more general) topological and functional charac- 
teristics may strongly alter the response of a network to 
thermal noise. Finding the general conditions required for 
reliable dynamics in a stochastic environment will be an 
important step towards a deeper understanding of the dy- 
namical features of real networks. 

We acknowledge the support of this work by the Hum- 
boldt Foundation. 
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